Automatic time stepping algorithms for implicit numerical simulations

An automatic time step size determination for non-linear problems, solved by implicit schemes, is presented. The time step calculation is based on the estimation of the integration error. This estimation is calculated from the acceleration difference. It is normalised according to the size of the problem and the integration parameters. This time step control algorithm modifies the time step size only if the problem has a long time physical change. On the other hand, Hessian matrix can be kept constant for several iterations however the problem is non-linear. According to the fact that the time step size is constant for some time step, the Hessian matrix shouldn’t be recalculated for each time step. A criterion selecting if Hessian matrix must be calculated or not is developed. Finally, a criterion of iterations divergence is also proposed. It avoids the determination, by the user, of a maximal iterations number. The iterations number is the smaller. Industrial numerical examples are presented that demonstrated the performances (precision and computational cost) of the algorithms. NOTATION. tn time after the time step number n h time step size q vector of nodal positions qn vector of nodal positions at time tn after convergence of the iterations qin vector of nodal positions at time tn after the iteration number i qn 0 vector of nodal positions at time tn after prediction qo vector of nodal positions at initial time. M mass matrix Fint vector of internal forces Fext vector of external forces ∃ first Newmark parameter ( second Newmark parameter ∀M first free parameter balancing sampling time around [tn, tn+1] for averaging inertia terms ∀F second free parameter balancing sampling time around [tn, tn+1] for averaging forces ∗ accuracy tolerance on the non-dimensional residue C tangent damping matrix K tangent stiffness matrix S Hessian matrix R residue of the equation of motion r non-dimensional residue of the equation of motion Τ pulsation for a one degree of freedom linear system Σ non-dimensional pulsation (Τh) ε integration error for a one degree of freedom linear system e integration error for the system PRCU tolerance of the integration error RAT factor multiplying the time step size between two steps INTRODUCTION. Non-linear dynamic problems can be solved with two kind of time-stepping algorithms: explicit or implicit. For an explicit code, solution at time tn depends only on solution at time tn-1. For an implicit code, solution at time tn depends on solution at time tn-1 but also at time tn. System is then solved with iterations. On the other hand, for an implicit algorithm the time step size can be longer than for an explicit algorithm. However the cost for a time step is higher for an implicit scheme, the time step number is smaller and the total time of calculation is often more interesting than for an explicit scheme. If the time step size is chosen too small, the calculation cost is very expensive. If it is chosen too big, the integration is not accurate enough or the iterations diverge. Therefor, time step size must be correctly calculated. According to the fact that the problem evolves with time, time step size must be adapted with this evolution. An automatic time step determination is then the only solution to solve accurately the problem in a short calculation time. For an industrial problem that has a large number of degrees of freedom, the most expensive operation of an implicit code is the inversion of the Hessian Matrix. For nonlinear problems, the Hessian matrix changes at each iteration. The Newton-Raphson iterations can sometimes converge when using the old inverted matrix. But this inverted matrix must be frequently recalculated otherwise the iterations diverge. This inversion occurs generally at begin of each time step and for some iteration selected by the user. But, if the Hessian matrix is not regularly inverted, problem diverges and if the inversion occurs too frequently, the problem becomes too expensive. According too the evolution of the problem with time, an automatic determination selecting if the inverted Hessian matrix must be recalculated or not, can reduce the total calculation cost. When the Newton-Raphson iterations diverge, the time step is rejected and the time step size is reduced. A problem is to determinate when the iterations diverge. Usually a maximum number of iteration is defined. If this number is too small a time step can be rejected however the problem slowly converges. If this number is chosen too big, some iteration are needlessly calculated when the divergence occurs. It is then interesting to determinate if divergence occurs in accordance with the residue evolution. The maximum number of iteration is more difficult to be correctly determinate when the inverted matrix is not calculated at each iteration. Indeed, this number depends on how frequently the inverted matrix is calculated. This paper proposes an automatic time step control algorithm based on the measure of the integration error. This algorithm modifies the time step size only if durable physical change occurs in the problem evolution. An automatic determination choosing if the Hessian matrix is recalculated is also proposed. This determination is based on residue evolution with iterations. Finally, a divergence criterion based on this residue evolution is implemented. Industrial numerical examples are then calculated with these new algorithms. NUMERICAL INTEGRATION OF TRANSIANT PROBLEMS WITH AN IMPLICIT CODE. The equations of motion of a non-linear structure, discretized in finite elements, can be write: ( ) ( ) 0 , , = − + q q F q q F q M ext int & & & & (1) First term of equation (1) is the vector of inertial terms. Second one is the vector of internal forces. This term is non-linear according to the geometrical non-linearity and plastic deformations. Last term is the vector of external forces. The non-linearity can be due from contact modelling, imposed displacements... The most general scheme for implicit integration, is based on averaged accelerations between time tn and tn+1. If positions and velocities at time tn are knew, positions and velocities at times tn+1 can be written: